Earth surface changes can be determined by comparing pairs of optical satellite images acquired on different dates. Precise images co-registration is a prerequisite in such applications, and this critical step is often a major source of limitation [1], [2]. For instance, a registration accuracy of less than ⅕ of a pixel is required to achieve a change detection error of less than 10% in Landsat Thematic Mapper images [3].
As to the measurement of ground surface displacements, most applications require a measurement accuracy of less than 1 m. This implies that the images co-registration accuracy should be even less, i.e., significantly smaller than the pixel size of most currently available optical satellite images. Examples of such applications include the measurement of coseismic ground deformations [4]-[7], ice flow [8], and sand dune migrations [9]. Difficulties in accurately co-registering satellite images arise from the non-ideal characteristics of the optical systems, the changing attitude of the spacecraft during the scanning operation of the images, digital elevation model (DEM) errors, and inaccurate resampling. The accuracy of the measurements of ground displacements, in addition, depends on the performance of the correlation technique. Despite these difficulties, encouraging results were obtained in a number of studies. It should be noted, however, that they were all carried on using data from only one imaging system and under restrictive conditions such as similar viewing angles and satellite tracks [4], [10], [11] or using external information from global positioning system (GPS) measurements [6]. Precise co-registration of images with viewing angle differing by more than 3° also seems out of reach [4], [11]. The operational use of such a technique, in particular to monitor coseismic deformations, would benefit from a more generic approach, allowing to cross-correlate images from different imaging systems with different viewing angles, and without the need for information other than what is extracted from the satellite ancillary data and the topography.
Known orthorectification and phase correlation models will be now briefly discussed.
A) Direct Orthorectification Model
The direct orthorectification model computes the geographic location on the ground where each pixel in the raw image, i.e., the focal plane of the instrument, has to be projected. Notations are derived from the SPOT satellite geometry handbook [15].
A1) Navigation Reference Coordinate System and Look Directions
The navigation reference coordinate system (O1,X1,Y1,Z1) is the spacecraft body fixed reference system.
Pushbroom satellite sensors consist of a CCD line array responsible for the image scanning operation. Expressed in the navigation reference coordinate system, the look directions are modeling the equivalent pointing direction of each CCD element. By being constant during the image acquisition, they provide the internal camera model accounting for the mirror rotation, optical distortions, and calibration parameters resulting from on-ground post-processing. The look directions are provided in ancillary data in the form of a two angle rotation (Ψx, Ψy) around the satellite body fixed system axes. See FIG. 1. For all columns c and for all rows r in the raw image, the look directions {right arrow over (u)}1 are given by
                                          u            →                    1                ⁡                  (                      c            ,            r                    )                    =                                                  u              →                        1            ′                    ⁡                      (                          c              ,              r                        )                                                                                                            u                  →                                1                ′                            ⁡                              (                                  c                  ,                  r                                )                                                          2                      ,          for      ⁢                          ⁢      all      ⁢                          ⁢      c        ,          r      =      1        ,    …    ⁢                  ,    N    with                                          u            →                    1          ′                ⁡                  (                      c            ,            r                    )                    =              (                                                                              -                  tan                                ⁢                                                                  ⁢                                                      Ψ                    y                                    ⁡                                      (                    c                    )                                                                                                                          tan                ⁢                                                                  ⁢                                                      Ψ                    x                                    ⁡                                      (                    c                    )                                                                                                                          -                1                                                    )              ,          for      ⁢                          ⁢      all      ⁢                          ⁢      r        ,  where N is the number of CCD elements in the line array.A2) Orbital Coordinate System and Attitude Variations
The orbital coordinate system (O2,X2,Y2,Z2) is centered on the satellite (O2=O1), and its orientation is based on the spacecraft position in space. See FIG. 2. Roll, pitch, and yaw variations are given as rotation angles around the Y2, X2, and Z2 axes defined by
         {                                                                                        Z                  →                                2                            ⁡                              (                t                )                                      =                                                            P                  →                                ⁡                                  (                  t                  )                                                                                                                                          P                      →                                        ⁡                                          (                      t                      )                                                                                        2                                                                                                                                      X                  →                                2                            ⁡                              (                t                )                                      =                                                                                V                    →                                    ⁡                                      (                    t                    )                                                  ⋀                                                                            Z                      →                                        2                                    ⁡                                      (                    t                    )                                                                                                                                                                                    V                        →                                            ⁡                                              (                        t                        )                                                              ⋀                                                                                            Z                          →                                                2                                            ⁡                                              (                        t                        )                                                                                                              2                                                                                                                                                          Y                    →                                    2                                ⁡                                  (                  t                  )                                            =                                                                                          Z                      →                                        2                                    ⁡                                      (                    t                    )                                                  ⋀                                                                            X                      →                                        2                                    ⁡                                      (                    t                    )                                                                        ,                              where {right arrow over (P)}(t) and {right arrow over (V)}(t) are the instantaneous position and velocity of the satellite, respectively, as shown in FIG. 2. The satellite look directions {right arrow over (u)}2(c,r) in the orbital coordinate system for all CCD elements are given, for all c, r=1, . . . , N, by {right arrow over (u)}2(c,r)=Rp(r)·Rr(r)·Ry(r)·{right arrow over (u)}1(c), where Rr(r), Rp(r), and Ry(r) are the roll, pitch, and yaw rotation matrices at the time of acquisition of image row r.A3) Look Directions in Terrestrial Coordinate System
For each pixel in the raw image, the corresponding look direction {right arrow over (u)}3 expressed within the terrestrial coordinate system is given by
                    u        →            3        ⁡          (              c        ,        r            )        =            [                                                                  X                                  2                  x                                            ⁡                              (                r                )                                                                                        Y                                  2                  x                                            ⁡                              (                r                )                                                                                        Z                                  2                  x                                            ⁡                              (                r                )                                                                                                        X                                  2                  y                                            ⁡                              (                r                )                                                                                        Y                                  2                  y                                            ⁡                              (                r                )                                                                                        Z                                  2                  y                                            ⁡                              (                r                )                                                                                                        X                                  2                  z                                            ⁡                              (                r                )                                                                                        Y                                  2                  z                                            ⁡                              (                r                )                                                                                        Z                                  2                  z                                            ⁡                              (                r                )                                                        ]        ·                                        u            →                    2                ⁡                  (                      c            ,            r                    )                    .      A4) Location on Earth Model
The corresponding ground location M of the raw image pixel (c,r) is determined by calculating the intersection between {right arrow over (u)}3(c,r) and the Earth ellipsoid model. For any of such pixel, the point M(xM,yM,zM) has to be found that verifies
                                                        O              3                        ⁢            M                    →                ⁢                  (                      c            ,            r                    )                    =                                                                  O                3                            ⁢              P                        →                    ⁢                      (            r            )                          +                  μ          ·                                                    u                →                            3                        ⁡                          (                              c                ,                r                            )                                            ,                  for        ⁢                                  ⁢        μ            >      0                                    and          ⁢                                          ⁢                                                    x                2                            +                              y                2                                                    A              2                                      +                              z            2                                B            2                              =      1        ,          with      ⁢                          ⁢              {                                                            A                =                                  a                  +                  h                                                                                                                          B                  =                                      b                    +                    h                                                  ,                                                        where O3 is the Earth Cartesian center and a and b are, respectively, the semimajor and semiminor axis of the ellipsoid. h is the approximated elevation above the ellipsoid at the ground location M.
Using a DEM, the intersection with the topographic surface is computed by locally and successively approximating the topography with a wider ellipsoid.
B) Phase Correlation Methods
Phase correlation methods have already shown good results for the measurement of ground deformation [4], [6], [7], [10]. All phase correlation methods rely on the Fourier shift theorem [23]: : The relative displacement between a pair of similar images is retrieved from the phase difference of their Fourier transform. Let i1 and i2 be two images that differ only by a displacement (Δx,Δy) such that i2(x,y)=i1(x−Δx,y−Δy).
By denoting by I1 and I2 their Fourier transform, from the Fourier shift theorem, we have the relationI2(ωx,ωy)=I1(ωx,ωy)e−j(ωxΔx+ωyΔy),where ωx and ωy are the frequency variables in column and row. The normalized cross-spectrum of the images i1 and i2 is then
                    C                              i            1                    ⁢                      i            2                              ⁡              (                              ω            x                    ,                      ω            y                          )              =                                                      I              1                        ⁡                          (                                                ω                  x                                ,                                  ω                  y                                            )                                ⁢                                    I              2              *                        ⁡                          (                                                ω                  x                                ,                                  ω                  y                                            )                                                                                                    I                1                            ⁡                              (                                                      ω                    x                                    ,                                      ω                    y                                                  )                                      ⁢                                          I                2                *                            ⁡                              (                                                      ω                    x                                    ,                                      ω                    y                                                  )                                                                  =              ⅇ                  j          ⁡                      (                                                            ω                  x                                ⁢                                  Δ                  x                                            +                                                ω                  y                                ⁢                                  Δ                  y                                                      )                                ,where * denotes the complex conjugate. The images' relative displacement can thus be estimated from the 2-D slope of the cross-spectrum's phase. By applying the inverse Fourier transform, we have the correlation functionF−1{ej(ωxΔx+ωyΔy)}=δ(x+Δx,y+Δy).
The images' relative displacement can then alternatively be estimated from the coordinates of the correlation peak. In case of subpixel displacements, this peak is not a Dirac delta function anymore, but a down-sampled version of a Dirichlet kernel [26]. Further processing is then required to recover the image shift. These approaches show that phase correlation methods fall into two categories. In the first category, the relative images' shift is recovered by explicitly estimating the linear phase of the images' cross-spectrum [4], [27], [28].
In the second category, the relative displacement is calculated by determining the exact location of the correlation peak [26]. This is generally not the case when images have been resampled for orthorectification. Also, to avoid correlation bias, frequency masking should be applied to only select parts of the cross-spectrum where the phase information is valid (images may be corrupted by aliasing or optical aberrations). For these reasons, a correlation algorithm whose main scheme belongs to the first category will be described, adaptive masking being applied on the cross-spectrum.
In [27], a robust approach has been proposed to evaluate the images phase difference. The normalized cross-spectrum matrix C(ωx,ωy) is, theoretically, a rank one matrix since C is separable. The idea of the study in [27] is to determine the best rank one approximation to the normalized cross-spectrum matrix. The displacement vector is recovered by calculating the slope of the unwrapped phase of the first singular vector, in each dimension. This method has proven a strong robustness against noise. However, there are two main drawbacks remaining. First, it is also subjected to phase wrapping. Even though this approach involves only 1-D unwrapping, it still remains a sensitive step. The second drawback, which is the main concern, is that the whole normalized cross-spectrum matrix (or a rectangular subset of it) has to be used to compute the best rank one approximation. This computation is potentially biased with corrupted phase values. A solution would be to use a weighted SVD, but most of these algorithms require the weight matrix to be positive definite symmetric [34]. Frequency weights with no a priori constraint on the spectrum orientation or separability should be applied.
In [4], another approach is proposed based on the Hermitian inner product of two functions. Define the theoretical normalized cross-spectrum of the images by C(ωx,ωy)=ej(ωxΔx+ωyΔy) and the one actually computed by Q(ωx,ωy). The projection of Q onto the continuous space defined by the theoretical cross-spectrums is defined as
                                          P                          Q              ,              C                                ⁡                      (                                          Δ                x                            ,                              Δ                y                                      )                          =                ⁢                              ∑                          ω              x                                ⁢                                    ∑                              ω                y                                      ⁢                                          Q                ⁡                                  (                                                            ω                      x                                        ⁢                                          ω                      y                                                        )                                            ⁢                                                C                  *                                ⁡                                  (                                                            ω                      x                                        ,                                          ω                      y                                                        )                                                                                            =                ⁢                              ∑                          ω              x                                ⁢                                    ∑                              ω                y                                      ⁢                                          Q                ⁡                                  (                                                            ω                      x                                        ⁢                                          ω                      y                                                        )                                            ⁢                                                ⅇ                                      -                                          j                      ⁡                                              (                                                                                                            ω                              x                                                        ⁢                                                          Δ                              x                                                                                +                                                                                    ω                              y                                                        ⁢                                                          Δ                              y                                                                                                      )                                                                                            .                                                        
The values of Δx and Δy that maximize the norm of this projection are the ones that are the most likely to solve the registration problem. It is then proposed to find (Δx,Δy) that maximizes the modulus |MPQ,C(Δx,Δy)|, where
                    MP                  Q          ,          C                    ⁡              (                              Δ            x                    ,                      Δ            y                          )              =                  ∑                  ω          x                    ⁢                        ∑                      ω            y                          ⁢                              M            ⁡                          (                                                ω                  x                                ,                                  ω                  y                                            )                                ⁢                      Q            ⁡                          (                                                ω                  x                                ,                                  ω                  y                                            )                                ⁢                      ⅇ                          -                              j                ⁡                                  (                                                                                    ω                        x                                            ⁢                                              Δ                        x                                                              +                                                                  ω                        y                                            ⁢                                              Δ                        y                                                                              )                                                                          ,and M(ωx,ωy) is a binary mask to filter out some frequencies. This technique is effective and insensitive to phase wrapping. It is suitable for both large and small displacement measurements; however, the resolution method proposed, based on a dichotomy, is computationally inefficient. Also, the frequency masking is not properly set.